Computing the current density \(j\) in the coil \(\def\curl{\operatorname{curl}}\def\Curl{\operatorname{Curl}}\def\div{\operatorname{div}}\)

Before we solve the non-linear magnetostatic problem, we first need to find the current \(j\) flowing inside the coil \(\Omega_c\).

The first property is that \(j\) is solenoidal, i.e. \(\div(j)=0\). Further, we presume the existence of a potential \(\phi\) with \(j=-\sigma\nabla\phi\), where \(\sigma\) denotes the connectivity. As for the boundary conditions, we prescribe \(n\cdot j = -\sigma\partial_n\phi=0\) on the exterior boundary \(\Gamma_{ex}\), and \(\phi=1\) and \(\phi=0\) on the inflow \(\Gamma_{in}\) and outflow boundary \(\Gamma_{out}\), respectively. Altogether, we have

\[\begin{split}\begin{align} j + \sigma\nabla\phi &= 0 \qquad\text{on }\Omega_c \\ \div j &=0 \qquad\text{on }\Omega_c \\ n\cdot j = \sigma\partial_n\phi &= 0 \qquad\text{on }\Gamma_{ex} \\ j &= 0 \qquad\text{on }\Gamma_{out} \\ j &= 1 \qquad\text{on }\Gamma_{in} \end{align}\end{split}\]

However, in our case, the coil is a loop with no inflow or outflow boundary. This is where the face “coil_cut_1” we defined in the geometry comes into play. For this purpose, we have to introduce “fictitious” points and introduce a clone of the face “coil_cut_1” in order to be able to prescribe the necessary boundary conditions.

We begin by loading the geometry and the mesh generated in the previous document

%%capture
%run TEAM_13_geometry.ipynb

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
File ~\AppData\Local\Temp\ipykernel_15972\1891279011.py:1
----> 1 import netgen.occ as occ
      2 from netgen.webgui import Draw as DrawGeo
      4 box1 = occ.Box(occ.Pnt(-100,-100,-50), occ.Pnt(100,100,50))

ModuleNotFoundError: No module named 'netgen'

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[3], line 1
----> 1 get_ipython().run_line_magic('run', 'TEAM_13_geometry.ipynb')

File ~\anaconda\envs\fem311\envs\mne\Lib\site-packages\IPython\core\interactiveshell.py:2456, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2454     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2455 with self.builtin_trap:
-> 2456     result = fn(*args, **kwargs)
   2458 # The code below prevents the output from being displayed
   2459 # when using magics with decorator @output_can_be_silenced
   2460 # when the last Python token in the expression is a ';'.
   2461 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~\anaconda\envs\fem311\envs\mne\Lib\site-packages\IPython\core\magics\execution.py:737, in ExecutionMagics.run(self, parameter_s, runner, file_finder)
    735     with preserve_keys(self.shell.user_ns, '__file__'):
    736         self.shell.user_ns['__file__'] = filename
--> 737         self.shell.safe_execfile_ipy(filename, raise_exceptions=True)
    738     return
    740 # Control the response to exit() calls made by the script being run

File ~\anaconda\envs\fem311\envs\mne\Lib\site-packages\IPython\core\interactiveshell.py:2981, in InteractiveShell.safe_execfile_ipy(self, fname, shell_futures, raise_exceptions)
   2979 result = self.run_cell(cell, silent=True, shell_futures=shell_futures)
   2980 if raise_exceptions:
-> 2981     result.raise_error()
   2982 elif not result.success:
   2983     break

File ~\anaconda\envs\fem311\envs\mne\Lib\site-packages\IPython\core\interactiveshell.py:294, in ExecutionResult.raise_error(self)
    292     raise self.error_before_exec
    293 if self.error_in_exec is not None:
--> 294     raise self.error_in_exec

    [... skipping hidden 1 frame]

File ~\AppData\Local\Temp\ipykernel_15972\1891279011.py:1
----> 1 import netgen.occ as occ
      2 from netgen.webgui import Draw as DrawGeo
      4 box1 = occ.Box(occ.Pnt(-100,-100,-50), occ.Pnt(100,100,50))

ModuleNotFoundError: No module named 'netgen'

Create the MESH object from the mesh generated by netgen

import sys
sys.path.insert(0,'../../../') # adds parent directory
import pde
MESH = pde.mesh3.netgen(geoOCCmesh)
MESH

np:4132, nt:23369, nf:3059, ne:720, nf_all:46966, ne_all:27728

The piece of code below duplicates the face, as described above, and generates a new MESH object.

import numpy as np

##############################################################################

face_index = pde.tools.getIndices(MESH.regions_2d,'coil_cut_1')
faces = MESH.f[MESH.BoundaryFaces_Region == face_index,:3]
new_faces = faces.copy()

points_to_duplicate = np.unique(faces.ravel())
new_points = np.arange(MESH.np,MESH.np+points_to_duplicate.size)

actual_points = MESH.p[points_to_duplicate,:]

t_new = MESH.t[:,:4].copy()
p_new = MESH.p.copy()
f_new = MESH.f.copy()

for i,pnt in enumerate(points_to_duplicate):

    # append point to list
    p_new = np.vstack([p_new,p_new[pnt,:]])

    # finding tets coordinates containing the ith point to duplicate
    tets_containing_points = np.argwhere(t_new[:,:4]==pnt)[:,0]

    for _,j in enumerate(tets_containing_points):
        #check if tet is left
        if MESH.mp_tet[j,0]<0:
            t_new[j,t_new[j,:]==pnt] = MESH.np + i

t_new = np.c_[t_new,MESH.t[:,4]]

for i,j in enumerate(faces.ravel()):
    new_faces.ravel()[i] = new_points[points_to_duplicate==j]

new_faces = np.c_[new_faces,np.tile(f_new[:,3].max()+1,(new_faces.shape[0],1))]
f_new = (np.r_[f_new,new_faces]).astype(int)

regions_2d_new = MESH.regions_2d.copy()
regions_2d_new.append('new')

identifications = (np.c_[points_to_duplicate,new_points]).astype(int)
# stop
MESH = pde.mesh3(p_new,MESH.e,f_new,t_new,MESH.regions_3d,regions_2d_new,MESH.regions_1d,identifications = identifications)
MESH

C:\Users\Radu\AppData\Local\Temp\ipykernel_8292\297165724.py:34: DeprecationWarning:

Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)

np:4149, nt:23369, nf:3076, ne:720, nf_all:47028, ne_all:27806

order = 1
D = pde.int.assemble3(MESH, order = order)
DB = pde.int.assembleB3(MESH, order = order)
N1,N2,N3 = pde.int.assembleN3(MESH, order = order)
unit_coil = pde.int.evaluate3(MESH, order = order, coeff = lambda x,y,z : 1+0*x, regions = 'coil')

face_in_1 = pde.int.evaluateB3(MESH, order = order, coeff = lambda x,y,z : 1+0*x, faces = 'coil_cut_1').diagonal()
face_in_2 = pde.int.evaluateB3(MESH, order = order, coeff = lambda x,y,z : 0+0*x, faces = 'coil_cut_1').diagonal()
face_in_3 = pde.int.evaluateB3(MESH, order = order, coeff = lambda x,y,z : 0+0*x, faces = 'coil_cut_1').diagonal()

###########################################################################

phi_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'M', order = order)
dphix_H1, dphiy_H1, dphiz_H1 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = order)
phiB_H1 = pde.h1.assembleB3(MESH, space = 'P1', matrix = 'M', shape = phi_H1.shape, order = order)

R0, RS0 = pde.h1.assembleR3(MESH, space = 'P1', faces = 'new,coil_cut_1')
R1, RS1 = pde.h1.assembleR3(MESH, space = 'P1', faces = 'coil_cut_1')

r = (face_in_1*N1 + face_in_2*N2 + face_in_3*N3) @ DB @ phiB_H1.T

M = phi_H1 @ D @ unit_coil @ phi_H1.T

K = dphix_H1 @ D @ unit_coil @ dphix_H1.T +\
    dphiy_H1 @ D @ unit_coil @ dphiy_H1.T +\
    dphiz_H1 @ D @ unit_coil @ dphiz_H1.T

r = -RS0 @ K @ R1.T @ (1+np.zeros(R1.shape[0]))
K = RS0 @ K @ RS0.T


RZ = pde.tools.removeZeros(K)
K = RZ @ K @ RZ.T

# M = RS0 @ M @ RS0.T
# M = RZ @ M @ RZ.T

r = RZ @ r

sigma = 1#58.7e6
from sksparse.cholmod import cholesky as chol
x = chol(sigma*K).solve_A(r)
x = RS0.T @ RZ.T @ x + R1.T @ (1+np.zeros(R1.shape[0]))

# MESH.pdesurf(x, faces = 'coil_face')

dx_x = (dphix_H1.T@x)*unit_coil.diagonal()
dy_x = (dphiy_H1.T@x)*unit_coil.diagonal()
dz_x = (dphiz_H1.T@x)*unit_coil.diagonal()

dphix_H1_P0, dphiy_H1_P0, dphiz_H1_P0 = pde.h1.assemble3(MESH, space = 'P1', matrix = 'K', order = 0)
unit_coil_P0 = pde.int.evaluate3(MESH, order = 0, coeff = lambda x,y,z : 1+0*x, regions = 'coil')
dx_x_P0 = (dphix_H1_P0.T@x)*unit_coil_P0.diagonal()
dy_x_P0 = (dphiy_H1_P0.T@x)*unit_coil_P0.diagonal()
dz_x_P0 = (dphiz_H1_P0.T@x)*unit_coil_P0.diagonal()

grid = pde.tools.vtklib.createVTK(MESH)
pde.tools.vtklib.add_H1_Scalar(grid, x, 'J')
pde.tools.vtklib.add_L2_Vector(grid,dx_x_P0,dy_x_P0,dz_x_P0,'grad_J')
pde.tools.vtklib.writeVTK(grid, 'current_density.vtu')

x.shape

(4149,)

import pyvista as pv
mesh = pv.read('current_density.vtu')
mesh

Header	Data Arrays
UnstructuredGrid Information N Cells 23369 N Points 4149 X Bounds -2.000e+02, 2.000e+02 Y Bounds -2.000e+02, 2.000e+02 Z Bounds -1.000e+02, 1.000e+02 N Arrays 3	Name Field Type N Comp Min Max J Points float32 1 0.000e+00 1.000e+00 Scalars_ Cells float64 1 0.000e+00 5.000e+00 grad_J Cells float32 3 -2.123e-03 2.110e-03

# mesh.plot(jupyter_backend='html')
# mesh.export_html('kek.html')
# p.add_mesh(mesh, style='wireframe', color='blue', label=None)
# p.add_mesh(mesh, label='Clipped')
# clipped = mesh.clip_scalar(scalars="Scalars_", value=1, invert=True)


mesh.set_active_scalars("Scalars_")
threshed = mesh.threshold([0,4])

# sample_function(
#     noise, [0, 1.0, -0, 1.0, 0, 1.0], dim=(20, 20, 20)

# clipped = mesh.threshold(value=1)
# # p.add_legend()
threshed.plot(jupyter_backend='html',show_edges=True,color="w")


# plotter = pv.Plotter()
# plotter.add_mesh(mesh.clip_scalar(scalars="Scalars_", value=100, invert=True), opacity=0.8)
# # plotter.show_axes()
# plotter.show(jupyter_backend='html')

pwd

'C:\\Users\\Radu\\Documents\\GitHub\\fem\\docs\\source\\notebooks'

import pyvista as pv
pv.set_jupyter_backend('trame')

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[4], line 2
      1 import pyvista as pv
----> 2 pv.set_jupyter_backend('trame')

File ~\anaconda\envs\fem311\Lib\site-packages\pyvista\jupyter\__init__.py:141, in set_jupyter_backend(backend)
     69 def set_jupyter_backend(backend):
     70     """Set the plotting backend for a jupyter notebook.
     71
     72     Parameters
   (...)
    139
    140     """
--> 141     pyvista.global_theme._jupyter_backend = _validate_jupyter_backend(backend)

File ~\anaconda\envs\fem311\Lib\site-packages\pyvista\jupyter\__init__.py:35, in _validate_jupyter_backend(backend)
     33 if backend not in ALLOWED_BACKENDS:
     34     backend_list_str = ', '.join([f'"{item}"' for item in ALLOWED_BACKENDS])
---> 35     raise ValueError(
     36         f'Invalid Jupyter notebook plotting backend "{backend}".\n'
     37         f'Use one of the following:\n{backend_list_str}'
     38     )
     40 # verify required packages are installed
     41 if backend == 'pythreejs':

ValueError: Invalid Jupyter notebook plotting backend "trame".
Use one of the following:
"ipyvtklink", "panel", "ipygany", "static", "pythreejs", "none"

sphere = pv.Sphere()
sphere.plot()